new method for numerical inversion of the Laplace transform. 



Bruno Hiipper and Eli Pollak 

Chemical Physics Department, Weizmann Institute of Science, 76100 Rehovot, Israel 

Abstract 

A formula of Doetsch (Math. Zeitschr. 42, 263 (1937)) is generalized and 
used to numerically invert the one-sided Laplace transform C((3). The neces- 
sary input is only the values of C((5) on the positive real axis. The method is 
applicable provided that the functions C((3) belong to the function space I? a 
defined by the condition that G(x) = e xa C(e x ), a > has to be square inte- 
grable. This space includes sums of exponential decays C(f5) = a n e~^ En , 
e.g. partition functions with a n = 1. 

In practice, the inversion algorithm consists of two subsequent fast Fourier 
transforms. High accuracy inverted data can be obtained, provided that the 
signal is also highly accurate. The method is demonstrated for a harmonic 
partition function and resonant transmission through a barrier. We find ac- 
curately inverted functions even in the presence of noise. 



I. INTRODUCTION 



It is often relatively easy to compute the Laplace transform 



C{(3) = / e~^ E C{E)dE 



(1.1) 



of a function rather than the function itself. Similarly, it is often known how to compute a 
function on the imaginary axis and it is desirable to have a useful method for analytic con- 
tinuation of the function to real time. Perhaps the most notable example is the computation 



which is straightforward in imaginary time t = —i%(3. A 'good' Laplace inversion method- 
ology would solve both of these issues. The difficulty is that the inverse Laplace transform 
is known to be an ill-posed problem, since the addition of a small perturbation (for example 



Different numerical methods have been worked out to attempt at overcoming this prob- 
lem 0,0. They divide roughly into five classes: The Fourier methods ||] which discretize 
the Bromwich inversion formula || 



This requires knowledge of the function in the complex plane and so does not really solve 
the problem. 

The next two classes are based on the idea that the original function C{E) may be 
expanded into a basis set of functions whose transforms are known. To this category belong 
the linear least-squares fitting methods, where different basis sets are used, e.g. orthogonal 
polynomials |l|||, sums of exponentials |7]||, rational fractions or continued fractions 0, 
as well as others ||10|| . Nonlinear least-square fits are necessary if the signal is decomposed 
directly into a sum of exponentials with unknown decay rates E n and coefficients a n |TT| , [12| . 
With both methods, it is difficult to treat signals accurately which are of the form C{(3) = 



of the propagator < x\e ltH / h \x' > which is very difficult because of the sign problem but 



{(3 — 1 — ib) 1 ) to the image C{(3) leads to a non-vanishing contribution (i.e. exp{{l + ib)E}) 
even in the limit of a very small perturbation (large b) Jl| . 
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Y^=i a n e~ En ^ . The Laplace transform of a polynomial-type basis possesses singularities and 
is inadequate for a fit to exponentials. On the other hand, an exponential with a non-integer 
decay-rate cannot correctly be approximated by exponentials of the sort e~ nl3 . As a result of 
these difficulties, these methods are able to give at most five exponentials. For other signals, 
such as rational polynomials, they have proved to be very accurate. 

Another approach is the singular value decomposition method (SVD) |l3| , p^| which is 
based on the theory of inverse problems. This method transfers the inverse Laplace transform 
into a matrix equation and the problem is transformed into the inversion of a nearly singular 
matrix, an ill-posed problem as well ||15||. 



The fifth and most recent approach is the maximum entropy method pJTql- In this 
method the entropy of the spectrum (which means in this context the number of ways of re- 
producing the spectrum) subject to a certain constraint is maximized. This approach allows 
to incorporate prior knowledge about the solution. Maximum entropy inversion is a nonlin- 
ear method, this complicates the practical application. However, it has proved its usefulness 
in recent computations, see for example Refs. [p|-|T9|j. The last two methods, maximum 
entropy and SVD, have recently been compared in simulating the electronic absorption 
spectrum of a chromophore coupled to a condensed phase environment and it turned out 
that the maximum entropy method is just able to reproduce the classical, smooth Franck- 
Condon contribution to the spectrum whereas SVD is capable of resolving finer oscillatory 



details @ 



In this paper we will resurrect an old formula, derived by Paley and Wiener |Tl|] and 



by Doetsch p2|| , which is direct, there is no need to use a basis set and there is no need to 
solve a set of nonlinear equations. The Paley and Wiener form was rederived by Gardner 
et al. |23| 40 years ago and applied with limited success (due in part to computational 
limitations) to the analysis of multi-exponential decay curves. The old formulae were derived 
for functions C(/3) which are L 2 integrable and so are not directly useful, for example for 
partition functions. We will generalize them, so that the method includes all functions which 
are L 2 a integrable, that is that the function G(x) = e xa C(e x ), a > is L 2 integrable. We 



find for an exponential series that the quality of the inversion depends on the magnitude of 
the n-th exponent E n . The smaller E n , the more accurate the inversion. This enables to 
enhance the resolution of the inverted data. 

In Section II, we derive the generalized Laplace inversion formula, numerical properties 
of the formula are discussed in Section III. The effect of shifting the signal is studied in 
Section IV. Applications to the harmonic oscillator partition function and a model resonant 
transmission probability are given in Section V. We end with a discussion of the merits of 
the method and outline some future extensions and applications. 

II. THE CONTINUOUS SET C' 1 OF INVERSE LAPLACE TRANSFORMS 

In this Section we derive and generalize a Laplace inversion formula which uses only the 
values of the Laplace transformed function C(/3) on the positive, real j3 axis. The starting 
point is the one-sided Laplace integral Eq. ( |1.1| ) for which we perform a transformation of 
variables 

E = S (3 = e x . (2.1) 

The motivation for this transformation, which goes back to Doetsch in 1936 [24|, is to 
substitute the Laplace kernel e~ l3E in which we have the product of the two variables by a 
different one which contains the sum of the two variables. As a result, the Laplace integral 
takes the form of a convolution. If on both sides of the Laplace integral transform the 
Fourier transform is applied, a certain convolution theorem can be used in order to express 
the right hand side of the integral equation as a product of two terms. Finally, an algebraic 
manipulation leads to the desired inversion formula. 

If we follow this route both sides of Eq. (|1 . 1| ) are multiplied with an exponential e xa 
with a > so that: 

/OO 
e a( x+ O e -e*+Z e f(l-a) C ( e {) ^ ( 2 .2) 
-oo L - 1 
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Now, the integrand on the right hand side consists of one part which depends only on the 
linear combination x + £ and a second, braced part which depends only on £. Next, both 
sides of the equation are Fourier transformed (with respect to x) and an application of the 
convolution theorem (which is equivalent to replacing the variable x by z = x + £) gives 



dxe lxy e xa C(e x ) = / 
The last integral can be written as 



3 C(i-«) C ( e f 



/ dz 



dze" e W 2y = / tit ^-i e -* = r(a + iy) 

-oo JO 



where r(x) denotes the Gamma function [p5 ]. Now, rearranging Eq. |2.2| leads to: 

1 



d£ 



3 «l-«) C ( e € 



dx e l ^e XQ C(e x ) 



r(a + iy) 

Fourier transformation of both sides of Eq. |2.5| yields the inversion formula 



gCC"- 1 ) pa g» 

C(£ = e 5 ) = o _ Jim / dy- 



2tt a^ooj_ a T(a + iy) 
Note that the inner integral in the inversion formula 



dx e lxy e xa C((3 = e x ) . 



9{y) 



dx < ir!l < '"(■{< ■'-. 



has the symmetry property 



g(-y) = g(y) 



and the Gamma function obeys 



t(z) = r z) 



(2.3) 



(2.4) 



(2.5) 



(2.6) 



(2.7) 



(2.8) 



(2.9) 



(where the bar denotes complex conjugation). This allows us to rewrite the inversion formula 
in a compact form as: 



C(E = e ? ) 



2tt 

7T 



lim / dy 



Re lim / dy 



v r(a + iy) 



g(y) 



;g{-y) 



o T(a + iy) J-oo 
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r(a - iy)" 
dx e ixy e xa C{e x ) 



(2.10) 



We have to require that e xa C(e x ) is square integrable, lest we encounter divergent inte- 
grals. This is not a very stringent requirement, as we can vary the parameter a to assure 
convergence. For example, the partition function of the harmonic oscillator 

m = 2^m (2 - n) 

leads to an L 2 integrable function provided that a > 1. 

Historically, Eq. |2.10| was first derived for a = 1 by Paley and Wiener |^TJ and then 
more rigorously by Doetsch p2| for the special case a = 1/2. In its form with a = 1 it 



was applied 40 years ago to the analysis of multicomponent exponential decay curves [23 



With the choice a > 1 the inversion formula is now amenable to a wider class of functions. 
Not less important is the fact that a proper choice of this exponent improves the numerical 
performance of the integrations by an order of magnitude as discussed in Appendix A. 

A generalization leading to the multi-dimensional inverse Laplace transform is straight- 
forward, as the main steps in the derivation of formula Eq. ( |2.10|) are based on the Fourier 
integral. The scalar variables E, f3 are replaced by iV-dimensional vectors E_, (3 to arrive at 
the iV-dimensional inversion formula: 

„£-(a-l) f a Jt-y r oo 

C{E) = 2 , . . Re lim / dy^ - -/ dx e^^C(/3) , (2.12) 

where the components of E_ are (e^ 1 , e^) and (3 = (e xi , e XN ). The components of a 
may be chosen to be different for each degree of freedom. 



III. NUMERICAL ANALYSIS 

In any numerical application two central questions arise: 

a) How does the accuracy with which C((3) is known, or which is reduced due to noise, affect 
the accuracy of the inversion technique? 

b) What is the range of (3 for which it is necessary to know the Laplace transformed function 
C(f3) in order to obtain a 'good' representation of the function C(E)7 
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To answer the first question we will consider in some detail the properties of the inversion 
formula for a single exponential C(/3) = e~ Eof3 . The original signal then is C(E) = 5(E — E ), 
where 5(x) is the Dirac '^'-function. The function g(y), cf. Eq. [2.7| , can be obtained 
analytically: 

g(y) = Eo a - iy na + iy). (3.1) 

This is a rapidly decaying function in the variable y, since the asymptotic behavior of the 
Gamma-function for fixed a and large \y\ is: 

\T(a + iy)\^ V2^\y\ a ~ 1/2 e- nW2 . (3.2) 

The asymptotic behavior of g(y) remains the same, even if the signal is a sum of many 
exponential terms. 

As seen from the exact integration Eq. ( ft.ip , the Gamma-function in the denominator of 
the inversion formula Eq. Q2.10| ) is cancelled and the y-dependence is simply an oscillating 
function with constant amplitude. This means that the numerically integrated function g(y) 
has to be exponentially small in y in order to compensate for the exponentially increasing 
factor coming from the denominator. The integrand of g(y) as shown in Fig. [I|, has an 
envelope of order unity and due to the term e lxy , it oscillates with increasing frequency as 
y increases. In the numerical integration, it is necessary to add up many Simpson boxes of 
height of the order of unity which leads to an exponentially small number e~ 7ry//2 . The boxes 
have different signs, so for large y this leads to the severe problem of adding many numbers 
with alternating signs, which differ only by a small amount. At this step we encounter the 
fact that the inverse Laplace transform is an ill-posed problem: a small error in the signal 

C{(3) =< C{(3) > +5C{(3) (3.3) 

leads to a large error in the integration 

< g(y) > +5g(y) = J K(x,y)(< C((3) > +6C((3)) . (3.4) 
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Let m denote the number of decimal digits of the signal and assume that the numerical 
integrations are all carried out with this same accuracy. Then g{y) can be obtained only up 
to a certain maximum value 

2 

2/max~-™hll0 (3.5) 
71 

beyond which any result is meaningless. This implies that the y-integration has to be 
truncated at y max - For the single exponent, this leads to the result: 

c(E) = 5m v;T o)> • (3 - 6) 

which is a function peaked at E — E with many side oscillations. An increase of the 
precision m allows to extend y max resulting in a narrower peak. In the limit y max — > oo the 
exact delta function is recovered. 

To appreciate this result it is instructive to compare it with an analogous analysis of the 
Fourier transform of a 5-function. In this case, truncation of the inverse Fourier transform 
leads to the approximation: 

C(E) = l si »(W£-£o)) (3J) 

7T E — Eq 

If the widths AE of these truncated 5-functions are defined as the distance between the first 
two zeros on each side of the maximum we have for the inverse Laplace transform 

7T 2,71 7T^ 

A L E = 2E sinh E ^0—^7: , (3.8) 

2/max ^rnax mill 10 

and for the inverse Fourier transform 

A F E = ^. (3.9) 

Whereas for the Fourier integral the width is independent of the position of the approximate 
delta, this is not true for the inverse Laplace transform. The resolution is automatically 
higher for lower energy. This means that for a given accuracy m the resolution will increase 
the lower the energy is. 
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Now, we come to the second task, the required range of j3. 

The integrand of g(y) is localized in x: for large negative x, the exponential factor e ax 
causes the integrand to decay, while for positive x, the exponentially decaying signal itself 
will also cause the integrand to decay rapidly. The decay for positive x is dependent on the 
magnitude of the lowest exponent E . If all calculations are performed with a precision of 
m digits, the value of the function is meaningful only up to the value /3 max defined as: 

e -£o/W = 1Q -m ^ /9Uc= ]n 10 (3.10) 

Eq 



Since the function g(y) is strongly decaying, the sampling theorem [^6| assures that 
the integrand needs to be sampled only at relatively few points. Therefore, the numerical 
integration involved in obtaining the function g(y) is straightforward. For example, for 
C{f3) = e~P, a = 1/2, the integrand is larger than 10 -15 in a range of Ax = 75.0. Using 
a Simpson integration rule with step size h = 0.1 leads to an integration error of ~ 10~ 16 
for y e [0,20]. That is, we have taken only 750 sampling points for this highly accurate 
integration, although the integrand varies quite rapidly for y = 20.0 with a period of A p x = 
0.3141. 

In summary, the range of j3 values needed as well as the resolution with which the inverse 
transform is performed is determined by the accuracy with which the Laplace transformed 
function C(f3) can be obtained. The lower the energy of the features, the larger the /3 range 
needed but also the higher the resolution. This property will be taken advantage of in the 
next section to significantly improve the algorithm. 

IV. SHIFTING THE LAPLACE TRANSFORM. 

Consider the case where the function C(E) is a sum of 5 functions: 

oo 

where the E/s are arranged according to ascending order. The Laplace transformed function 
is 
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(4.2) 



3=0 



Define the shifted Laplace transformed function as: 



C{(3,E s )=eP E °C{(3). 



(4.3) 



The inverse Laplace transform will lead to the function, 



oo 



C(E,E s ) = Y,a J S[E-(E 1 -E s )], 



(4.4) 



where all the eigenenergies have been moved closer to the origin by the amount E s . As 
shown in the previous section, such a shift will lead to enhanced resolution in the inverted 
function. 



we plot the inverse transform with a cutoff at y max = 5.0, which means that the accuracy 
of the signal is only 3 decimal digits. Even the lowest decay rate can hardly be estimated 
accurately as may be inferred more clearly from a blow-up of the dashed line shown in Fig. 
|3|. The width of the lowest 5-f unction is of the same order as the spacing and so it is hardly 
discernible. Shifting the function by E s = 0.9 gives a dramatic increase in resolution. A 
blow-up of this first peak is provided in Fig. f|. From this figure we can find that the 
maximum lies at E = 0.0995. The price to be paid for the increased resolution is that f3 max 
(cf. Eq. p,10| ) must be increased, since it is inversely proportional to the magnitude of the 
lowest eigenvalue which has now been reduced from Eq to Eq — E s . One may now repeat the 
computation, shifting the data by 0.999 instead of 0.9 and the peak will be resolved with 
even higher accuracy. In this way, the eigenvalue can be obtained with arbitrary accuracy. 



For illustration, let us consider four exponentials with decay rates 1,2,3 



,4. In Fig. | 



V. NUMERICAL APPLICATIONS 



A. Partition function of the harmonic oscillator 



The exact inversion of the partition function Eq. (|2.11 ) leads to a train of delta functions 
at the positions of the eigenvalues of the harmonic oscillator. This function was chosen 
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because its numerical Laplace inversion belongs to the most difficult class of problems. A 
non-linear least-squares method (without any knowledge in advance) could fit at most five 
exponentials. On the other hand, expansions in different basis sets converge too slowly 

0. The inverse Laplace transform of the partition function was computed with different 
degrees of decimal digits precision. Fig. [5] compares calculations with double precision, 

1. e. 15 decimal digits, and a little higher accuracy, 26 decimal digits. Whereas for double 
precision only the two lowest eigenvalues can be identified, at the higher accuracy the four 
lowest eigenvalues are resolved. 

The results of pushing the accuracy to 60 and 105 decimal digits precision are shown 
in Fig. H At 105 decimal digit precision it is possible to identify the eigenvalues up to 
the 10th level. The range of f3 values used in all these computations is as in Eq. fl3.10|) , 
Anax ~ 4.5m. Of course, these calculations cannot be applied to data obtained from a Monte 
Carlo computation. However, as also discussed in the next section, they may be used to 
invert basis sets which can then be fitted to Monte Carlo data. These results also serve to 
demonstrate the relative simplicity and accuracy of the method and the fact that in principle 
it will work for any number of peaks. 

To test the noise-sensitivity of the inverse Laplace transform, we added to the signal a 
Gaussian distributed noise with zero mean and different levels of RMS deviation a. The 
signal is assumed to be given up to x = 5.52, i.e. /9 max = 250. Fig. [7] shows that beyond the 
cut-off value 2/ max there is an accumulation of numerical errors and the signal deviates from 
a cosine-like wave. This Figure also confirms that the cut-off value depends rather linearly 
on the logarithm of the RMS deviation of the noise a. In Fig. §, the signal is shifted to 
the left by E s = 0.4, so that the smallest decay rate is around E m 0.1. The cut-off values 
change only slightly under the shift operation, but the integrand contains more oscillations 
before the cut-off, leading to an enhanced resolution in the peaks. 
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B. Reflection probabilities 



The Laplace transform of the reflection probability for the Eckart potential [28 



R(E) = 1 + oosh(V4a2 - ^ ^ 

cosh(2aJE/V$) + cosh( v / 4a 2 - tt 2 ) 



see Fig. is computed by numerical integration. Then the real inversion formula Eq. ( |2.10| ) 
is used to regain the reflection probability. The difference between the exact function and 
the inverted one for the parameter choice a = 4.0, V* = 5.0 is too small to be seen by the 



naked eye. A blow-up of the error is shown in Fig. |1C| Even for the rather low accuracy 
of only 3 decimal digits the relative error is about 10~ 2 , and as seen from the Figure, it 
decreases with increasing precision of the data. For the parameter a = 12.0, = 5.0 the 



results are a bit worse, as shown in Fig. [IT], due to the 'Gibb's phenomenon' ||29|| . Near the 
step, E pa 5.0 the error increases significantly. 

In all the computations the cut-off y max was chosen to minimize the error: decreasing the 
value of y max reduces the resolution but increasing it leads to numerically wrong values due to 



the uncertainty of the signal. In Fig. [12] we show a typical integrand Ke{g(y)/T (1/2 + iy)}: 
if the inverse Laplace transform is not known, it is easy to judge which value for y max has 
to be chosen, as the integrand decays smoothly and then produces artificial oscillations and 
blows up. (For an exact step transmission probability f((3) = ^e~ vt/3 , the integrand goes 
asymptotically to as l/(a — 1 + iy).) 

C. Below-barrier resonance 

A small resonance in the form of a Lorentzian is added to the transmission probability 

T(E) g2 | coshE/lW 

1 ; [E-E ) 2 + e 2 coshE/Vo + a' { ' } 



with parameters given in Fig. [U| The accuracy of the data is 10~ 6 and the features 
are reproduced quite well. The oscillations at very low energy are side oscillations of the 
resonance and it is possible to smooth them away. As outlined above, the resolution depends 
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on the energy of the feature. In order to reproduce a Lorentzian of width e, it is necessary 
to have at least a comparable resolution. To check whether the Lorentzian coming out of 
the inversion is broadened because of lack of resolution, the signal can be shifted towards 
lower energy. 

In this example we took the Laplace transform of the reflection probability with /5 max = 
10 6 . One may also use the transmission probability, however it diverges at (5 — and so this 
requires some care. 

VI. DISCUSSION 

In this paper we have resurrected and generalized a formula of Doetsch which enables a 
direct Laplace inversion of a large class of functions. By suitable scaling, these can include 
functions that are not L 2 integrable. Therefore the algorithm is directly applicable to parti- 
tion functions, for example. The method is relatively simple, all that is needed are two fast 
Fourier transforms. It is not necessary to pre-smooth the data. The method is controllable, 
the more accurate the Laplace inverted data, and the larger the range, the more accurate 
are the inversion results. The parameters of the inversion are controlled by the accuracy of 
the data only. As a result, the method is stable with respect to small perturbations. 

We have shown that in practice, an extremely high quality inversion can be obtained 
provided that the signal is also of very high accuracy. This is not merely an academic 
exercise. For example, the Laguerre basis set may be taken, systematically inverted, and 
the resulting numerical functions may be stored. Then the Laplace transformed function 
may be expanded in terms of Laguerre polynomials. The inverted function is then obtained 
merely by reading off the inverted Laguerre functions. The utility of such a procedure 
depends on the qualitfy of the fit of the polynomials to the numerical Laplace transformed 
data. It may be, that more sophisticated techniques should be used which include local 
smoothing of the data, such as the DAFS methodology |[33|| . In any case, once the Laplace 
transformed data is projected onto standard basis sets, the high accuracy inversion may be 
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used to obtain the inverted function. 

An important property of the inversion technique is the fact that the resolution of the 
resulting signal depends on the location of the signal. The closer it is to the origin, the higher 
is the resolution. This allows for a shifting of the signal to obtain an increased resolution. 
The price to be paid is that each shift demands knowledge of the function for larger values 
of /3. For analytic functions, such as the Laguerre polynomials, this does not present any 
severe difficulty, as present day computers enable computations with very high accuracy, 
which is also demonstrated for the harmonic oscillator partition function. 

The Laplace inversion method presented in this paper is ideally suited for data obtained 
from matrix multiplication schemes f31| , |32|j . These methods produce the data at points 
Pj = ||30|| , while the inversion requires f3j = e^ Ax . 

In this paper we have not considered correlation functions. Elsewhere [[34] we will present 



the application of the present method to spectra and correlation functions. In principle there 
is no special complication except for the fact that in some cases a two dimensional inverse 
Laplace transform has to be computed. 

We have also not considered directly the numerical analytic continuation of functions. 
As already mentioned in the Introduction, once the inverted function is obtained, it may 
be Fourier transformed to obtain the analytically continued function. In this sense, the 
inversion technique presented in this paper may be thought of as a representation of the 
complex valued Dirac 5 function. The real question is one of practical usage, that is the 
level of accuracy needed to obtain the real time function from the imaginary time function 
as well as the range of (3 values needed for a given time length. Other applications are 
the computation of moments of a probability distribution from its transform ||35|| . These 
questions will be considered in future studies |34[ . 
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APPENDIX A: OPTIMIZING THE CHOICE OF a 

We will outline how the parameter a can help reduce the numerical effort drastically, 
especially in high precision calculations. 

The main numerical advantage of introducing a is a shortening of the integration interval 
in x needed for obtaining g(y), cf. Eq. (|2.7| ). The range of integration [x m i n , x max ] is 
determined by the required accuracy s = 10~ m . The negative limit is mainly fixed by the 
exponential e ax , 

1 e 

Xmm = — hi- , (6.1) 

a C(0) 

and the positive limit is due to the very rapid decay of C(e x ) which is almost independent 
of a and is determined by the smallest decay rate. The larger a, the smaller the integra- 
tion interval, but if a becomes too large the integrand increases exponentially, magnifying 
uncertainties in the signal. 

The maximum value of the integrand, if one exponential decay C(j3) = a§e~ E °^ is con- 
sidered, is at x m = \na/E and the integrand I(x) takes the value 

I{x m ) = a e a{hia/E °- 1) w a e alna , (6.2) 

which goes essentially as a!. The larger a, the more digits are required in the computation. 
On the other hand, the outcome of the integration must cancel the denominator r(a + iy) 



whose large y-asymptotics is given by Eq. ( |3.2|) . For a given y max the order of magnitude of 



r(a + iy) divided by the integrand at y = 0, T(a) ~ a a rs e alna h as to be comparable to 
the given accuracy e = 10 _m : 

a-l/2 -7ry max /2 p 

1 x\\ = — 6 - 3 

(a — 1)1 a 

In summary, for large cut-off values y max the stepsize in the x integration remains ap- 
proximately the same. A change of a reduces the interval of the first integration, but to 
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keep the same resolution (i.e. keep y max fixed) it is necessary to increase the precision m. 
We found that for m 100, a 10 is a reasonable choice. 
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FIGURES 



FIG. 1. Integrand of the inverse Laplace inversion formula for a signal of one exponential decay. 
The envelope is decaying exponentially for x —* — oo and even more rapidly for x — > +00. The 
rapid oscillations result in an exponentially small value g(20) ~ 5 • 10~ 14 of the integral although 
the integrand is of the order of unity. 

FIG. 2. Inverse Laplace transform of a sum of four exponential decays with decay rates 
E n = 1, 2, 3, 4. The accuracy of the signal is taken as 3 decimal digits. The inversion of the original 
data allows at most the estimation of the first delta function (dashed line) at E = 1. The solid 
line shows the inversion of the data shifted by 5E = 0.9 to the left. The first maximum can now 
be estimated much more accurately. 

FIG. 3. Magnification of the unshifted inversion of Fig. The exact curve should yield 
a delta function at E = 1.0. Due to the insufficient accuracy of the data, the four components 
overlap and distort the maximum to E ~ 1.05. 

FIG. 4. Magnification of the shifted inversion Fig. |^. The exact curve should yield a delta 
function at E = 0.1. This value may now be estimated very accurately from the shifted data, even 
though the accuracy (m = 3) is low. 

FIG. 5. Numerical inverse Laplace transform for the partition function of the harmonic oscil- 
lator. The exact inverse should yield J2n^(^ — (n + 1/2)). The two lines correspond to different 
input signals whose accuracy (significant decimal digits) is indicated in the insert. The value of 
a = 4 was used for all computations with the harmonic oscillator partition function. 

FIG. 6. High precision numerical inverse Laplace transform for the partition function of the 
harmonic oscillator. Other notation is as in Fig. ||. 

FIG. 7. Noisy data. The integrand of the real inversion formula for the partition function 
of the harmonic oscillator is plotted vs. y. Gaussian noise with RMS deviation a as indicated is 
added to the signal and this leads to a reduction of the cut-off value for the y-integration. 

FIG. 8. Noisy shifted data. The data used for Fig. [7] are shifted by E s = 0.4 to the left. 
The cut-off values in y remain the same, but because of the faster oscillation of the integrand, the 
resolution of the final inversion peaks will be increased. 

FIG. 9. Reflection probabilities for the Eckart barrier with two different choices of the param- 
eters. For all reflection probabilities we used a = 0.5. 
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FIG. 10. Logarithm of the error of the inverted reflection probability of the Eckart potential 
with a = 4.0, = 5.0. The signal for the inversion is obtained by numerical Laplace transform of 
the exact reflection probability and the accuracy in decimal digits of the numerical Laplace integral 
is indicated. The values for y max are 5.5 and 12.0 for 3 and 8 digits accuracy respectively. 



FIG. 11. Logarithm of the error of the inverted reflection probability of the Eckart potential 
with a = 12.0, V* = 5.0. Other notation is as in Fig. 
E = 5.0 due to Gibb's phenomenon. 



10. The error increases near the step at 



FIG. 12. Integrand of the real inversion formula for the Eckart barrier reflection probability 
at a = 4.0, V* = 5.0. The integrand is expected to decrease like l/(c + iy), but beyond the cut-off 
2/max ~ 12.5 artificial oscillations arise and the integrand blows up. 



FIG. 13. Numerical inverse Laplace transform for a below barrier resonance added to the 
reflection probability of the Eckart barrier T res (E) = e 2 /((E-E ) 2 +s 2 ), with e = 0.013, E = 0.05, 
added to the transmission probability T(E) = (cosh20i? — 1)/(100 + cosh20-E). The accuracy of 
the data is 6 decimal digits. 



20 




21 




22 




Fig. 4 

23 




24 




Fig. 6 



25 




26 



+ 



250.0 



150.0 




-150.0 



8.0 10.0 



Fig. 8 



27 




20.0 



E 



Fig. 9 



28 




29 



I 



o 




10.0 



Fig. 11 



30 



5.0 




Fig. 12 



31 




32 



